#include "../DWMRICore/Fiber.h"
#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/Point3I.h"

#include <time.h>
#include <omp.h>

char *fiber_pathname = NULL;
char *dst_pathname = NULL;

char **confPathname = NULL;
char **vertPathname = NULL;
int w, h, d, scale;
float stepSize, maxLength;

float fiber_length = 0.0f;

void PrintUsage()
{
	printf("Usage: DWILIC -fiber <configure_file> -dst <dst_file> -length <fiber_length>\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-fiber") == 0) {
			fiber_pathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			dst_pathname = argv[++i];
		} else if (strcmp(argv[i], "-length") == 0) {
			fiber_length = (float)atof(argv[++i]);
		} else {
			return -1;
		}
	}

	if (fiber_pathname == NULL || dst_pathname == NULL || fiber_length < 0.0001f)
		return -2;

	return 0;
}

float FindNextVoxel(CPoint3F &pos, CVector3F &dir, CPoint3I voxel, CPoint3I &next)
{
	float u_x = (voxel.x + 0.5f) / scale;
	float u_y = (voxel.y + 0.5f) / scale;
	float u_z = (voxel.z + 0.5f) / scale;
	float l_x = (voxel.x - 0.5f) / scale;
	float l_y = (voxel.y - 0.5f) / scale;
	float l_z = (voxel.z - 0.5f) / scale;

	float min_l = 1000.0f;
	float l = 0.0f;

	// case 1
	l = (u_x - pos.m_x) / dir.m_x;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x + 1;
		next.y = voxel.y;
		next.z = voxel.z;
	}

	// case 2
	l = (l_x - pos.m_x) / dir.m_x;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x - 1;
		next.y = voxel.y;
		next.z = voxel.z;
	}

	// case 3
	l = (u_y - pos.m_y) / dir.m_y;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x;
		next.y = voxel.y + 1;
		next.z = voxel.z;
	}

	// case 4
	l = (l_y - pos.m_y) / dir.m_y;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x;
		next.y = voxel.y - 1;
		next.z = voxel.z;
	}

	// case 5
	l = (u_z - pos.m_z) / dir.m_z;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x;
		next.y = voxel.y;
		next.z = voxel.z + 1;
	}

	// case 6
	l = (l_z - pos.m_z) / dir.m_z;
	if (l > 0.0f && l < min_l) {
		min_l = l;
		next.x = voxel.x;
		next.y = voxel.y;
		next.z = voxel.z - 1;
	}

	min_l += 0.0001f;

	pos.m_x += dir.m_x * min_l;
	pos.m_y += dir.m_y * min_l;
	pos.m_z += dir.m_z * min_l;

	//voxel.x = next.x;
	//voxel.y = next.y;
	//voxel.z = next.z;

	return min_l;
}

float ComputeValue(CScalar *noise, CPoint3F v1, CPoint3F v2)
{
	CPoint3F start = v1;
	CPoint3F end = v2;
	CVector3F dir = end - start;

	float length = 0.0f;
	float totalLength = dir.Length();
	dir.Normalize();

	CPoint3I voxel;
	voxel.x = int(floor(v1.m_x * scale + 0.5f));
	voxel.y = int(floor(v1.m_y * scale + 0.5f));
	voxel.z = int(floor(v1.m_z * scale + 0.5f));
	
	float v = 0.0f;
	
	while(length < totalLength) {
		CPoint3I next;
		float thisLength = FindNextVoxel(start, dir, voxel, next);

		if (thisLength > totalLength - length) {
			v += (totalLength - length) * noise->GetScalarFast(voxel.x, voxel.y, voxel.z);
			break;
		} else {
			v += thisLength * noise->GetScalarFast(voxel.x, voxel.y, voxel.z);
		}

		length += thisLength;
		voxel = next;
	}
	
	return v;
}

int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	/* read the configure file */
	ReadConfigureFile(fiber_pathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);

	/* generate the white noise volume */
	clock_t time = clock();
	srand(time);
	CScalar *noise = new CScalar();
	noise->CreateScalar(w, h, d);
	for (int i = 0; i < w*h*d; ++i) {
		float temp = (float)rand();
		noise->SetScalar(i, temp);
	}

	float length = fiber_length > maxLength ? maxLength : fiber_length;
	uint count = uint(length / stepSize);
	int cw, ch;
	CFiber *fibers = NULL;

	CScalar *lic = new CScalar();
	lic->CreateScalar(w, h, d);

	omp_set_num_threads(64);

	for (int z = 100; z < 150; ++z) {
		ReadFibers(confPathname[z-100], vertPathname[z-100], &fibers, cw, ch);

		for (int y = 0; y < h; ++y) {
#pragma omp parallel for   
			for (int x = 0; x < w; ++x) {
				int index = y * w + x;


				if (fibers[index].m_fCount + fibers[index].m_bCount < 3) 
					continue;

				float v = 0.0f;
				uint fc = fibers[index].m_fCount > count ? count : fibers[index].m_fCount;
				uint bc = fibers[index].m_bCount > count ? count : fibers[index].m_bCount;
				for (uint k = 1; k < fc; ++k) {
					v += ComputeValue(noise, fibers[index].m_pF[k-1], fibers[index].m_pF[k]);
				}
				for (uint k = 1; k < bc; ++k) {
					v += ComputeValue(noise, fibers[index].m_pB[k-1], fibers[index].m_pB[k]);
				
				}

				lic->SetScalar(x, y, z, v);
			}
		}

		SafeDeleteArray(fibers);

		printf("z = %d\n", z);
	}

	lic->SaveScalarFile(dst_pathname);

	SafeDelete(noise);
	SafeDelete(lic);

	return 0;
}